      doubleprecision Function myrotlaw(r,tau,beta) result(chi)
C     r - cylinder radius, r < 1 inside the star, 
C     0<=tau<=1 - time of collapse, 
C     0<=beta<=1 - patameter of rotation
C     chi=rot. potential
      Implicit none
	integer, parameter :: Ndots=100
      real(8) r,tau,beta,solve_rot,rho_polytrope,find_average
	real(8) tau1,beta1,r_loc(Ndots),r_loc2(Ndots)
      real(8) ksi(Ndots),rho(Ndots)  ! "radius" and average density of initial star (assume polytrope n=3)
	real(8) pi,x0,y,ksil,yl,ksir,yr,ksi0,y0
	real(8) t1,t2,r2
	integer k,km1,kp1
	real(8), save :: Xdot(Ndots),Chidot(Ndots),tau_save,beta_save
	data pi/3.14159265358979323846264338328D0/
	data tau_save/10.D0/beta_save/10.D0/
c      data ksi/0.D0,0.05D0,0.1D0,0.16D0,0.25D0,0.35D0,0.5D0,0.61D0,
c     1     0.75D0,0.86D0,1.D0,1.11D0,1.25D0,1.37D0,1.5D0,1.75D0,2.D0,
c     2     2.25D0,2.5D0,2.75D0,3.D0,3.25D0,3.5D0,3.75D0,4.D0,4.25D0,
c     3     4.5D0,4.75D0,5.D0,5.25D0,5.5D0,5.75D0,6.D0,6.25D0,6.5D0,
c     4     6.7D0,6.8968475D0/
c      data rho/0.333333333333333D0,0.33308215220D0,0.3323357753D0,
c     1	 0.3307881371D0,0.3271706489D0,0.3214145343D0,0.309679049D0,
c     2	 0.2990361166D0,0.2835083979D0,0.2701320206D0,0.2521292237D0,
c     3     0.2375669019D0,0.2189323887D0,0.2031774259D0,0.1866077497D0,
c     4     0.1569132288D0,0.130745426D0,0.1083557687D0,8.95875502D-2,
c     5     7.406970909D-2,6.134994017D-2,5.097414763D-2,4.25273668D-2,
c     6     3.565033531D-2,3.004226858D-2,2.545639767D-2,2.169309757D-2,
c     7     1.859206739D-2,1.602520362D-2,1.38902945D-2,1.210581405D-2,
c     8     1.060663739D-2,9.340642717D-3,8.26598585D-3,7.348971975D-3,
c     9     6.710378376D-3,6.152050295D-3/

	if(tau.gt.1.D0.OR.tau.lt.0.D0.OR.beta.lt.0.D0.OR.beta.ge.1.D0)then
	 print*,'Wrong values in myrotlaw'
	 print*,'beta=',beta
	 print*,'tau=',tau
	 read*
       stop
      end if

c----/calculate the distribution of dots along 0<r<1/------
      call make_dots(Ndots,1.d-3,ksi)
    
c----/find average rho in polytrope n=3------     
	do k=1,Ndots
	rho(k)=rho_polytrope(ksi(k))
	end do

      beta1=beta*0.5D0*rho(Ndots)
      tau1=tau*pi/2.D0/dsqrt(rho(1)*(1-beta1/rho(1))**3)

      if(tau.ne.tau_save.OR.beta.ne.beta_save)then ! initialize for the first time usage
	 tau_save=tau
	 beta_save=beta
	 do k=Ndots,1,-1
	  if(k.ne.Ndots)then
	   x0=Xdot(k+1)
        else
	   x0=0.99d0
	  end if 
	  Xdot(k)=solve_rot(x0,tau1*dsqrt(rho(k)),beta1/rho(k))
	 end do
       
	 r_loc=ksi*Xdot
	 r_loc=r_loc/r_loc(Ndots)
	 r_loc2=r_loc**2

	 Chidot(1)=0.D0
	 do k=2,Ndots
	  km1=k-1
	  x0=find_average(Xdot(k),Xdot(km1))
	  Chidot(k)=Chidot(k-1)+x0*(r_loc2(k)-r_loc2(km1))
	 end do
	end if

      if(r.ge.1.D0)then
	 chi=Chidot(Ndots)
	 return
      elseif(r.eq.0.d0)then
	 chi=0.d0
	 return
	end if

      k=minloc(r-r_loc,mask=r-r_loc.ge.0.D0,dim=1)  
	kp1=k+1

	y=r*Xdot(Ndots)
	yl=Xdot(k)*ksi(k)
	yr=Xdot(kp1)*ksi(kp1)
      ksil=ksi(k)
	ksir=ksi(kp1)
	do
	 ksi0=(y-yl)*(ksir-ksil)/(yr-yl)+ksil
       t1=rho_polytrope(ksi0)
       x0=solve_rot(y/ksi0,tau1*dsqrt(t1),beta1/t1)
	 y0=x0*ksi0
	 if(dabs((y0-y)/y).le.1.d-9)then
	  exit
       elseif(y0.le.y)then
	  yl=y0
	  ksil=ksi0
       else
	  yr=y0
	  ksir=ksi0
	 end if
      end do
    
      chi=Chidot(k)
      r2=r**2
      t1=(r2-r_loc2(k))/(r_loc2(kp1)-r_loc2(k))

	t2=find_average(Xdot(k),x0)
      chi=chi+t1*t2*(r_loc2(kp1)-r2)

	t2=find_average(Xdot(k),Xdot(kp1))
      chi=chi+t1*t2*(r_loc2(kp1)-r_loc2(k))

	t2=find_average(x0,Xdot(kp1))
      chi=chi+t1*t2*(r2-r_loc2(kp1))

      return
	end function
C=====================================================================================================
C*****************************************************************************************************
      doubleprecision function solve_rot(x0,tau,beta) result(x)
C     solve the equation for x=r/r0 (the compression), using local tau and beta
C     x0 is an initial guess about x
      implicit none
	real(8) rho,tau,beta
	real(8) x0,x1,x2,xmin,y,dy

	xmin=beta/(1-beta)

	if(x0.le.xmin.OR.x0.ge.1.d0)then
	 x=0.5D0*(xmin+1)
      else
	 x=x0
      end if

	do 
	 x2=1-beta
	 x1=dsqrt((1-x)*(x*x2-beta))
       y=x1/x2+dasin(dsqrt(x2*(1-x)/(1-2*beta)))/x2**1.5D0

	 if(dabs(y-tau).le.1.d-9)exit
	 dy=-x/x1
	 x1=x+(tau-y)/dy
	 if(x1.ge.1.D0)then
	  x=0.5D0*(x+1)
       elseif(x1.le.xmin)then
	  x=0.5D0*(x+xmin)
       else
	  x=x1
       end if
	end do
    
      return
      end function

c***********************************************************************************************
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine make_dots(N,delta,dots)
C     to calculate the distribution of dots in the domain 0<rad<1
C     N is the number of dots
C     delta is the value of first step, x_{k}=x_{k-1}*alpha+delta,
C     alpha^N-1=(alpha-1)/delta
      implicit none
	integer N,i
	real(8), dimension(N) :: dots
	real(8) delta,alpha,alpha_old

	if(N*delta.ge.1.d0)then
	 print*,'Forbidden values in make_dots',N*delta
	 read*
	 stop
      end if

      alpha=1-(delta+1.d0/dble(N))*dlog(N*delta)
	alpha_old=alpha
      
      do
	 alpha=(1+(alpha-1)/delta)**(1.d0/dble(N))
	 if(dabs(alpha-alpha_old).le.1.d-10)exit
	 alpha_old=alpha
	end do

      dots(1)=0.d0
      do i=2,N
	 dots(i)=dots(i-1)*alpha+delta
	end do

	dots=dots/dots(N)
	return
	end subroutine

c***********************************************************************************************
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      doubleprecision function rho_polytrope(x) result(rho)
c     to calculate average rho=const*M/R^3 
c     for polytrope n=3 (gamma=4/3) (INDX=1) or
C     for polytrope n=3/2 (gamma=5/3) (INDX=2) or
C     for polytrope n=1 (gamma=2) (INDX=3)
c     as a function of dimensionless radius 0<x<1. rho(0)=1
c     the function rho(x) is approximate!
      implicit none
	integer, parameter :: INDX=2
	real(8) x,x2
	real(8) AV(4,2)
	data AV(:,1)/32.80686d0,922.98808d0,-342.53735d0,40567.28285d0/
	data AV(:,2)/4.899943d0,22.913686d0,-28.217443d0,112.4382012d0/

      x2=x**2
      if(INDX.eq.1.or.INDX.eq.2)then
	rho=1+x2*(AV(1,INDX)+x2*(AV(2,INDX)+
     *	    x2*(AV(3,INDX)+x2*AV(4,INDX))))
	rho=1.d0/rho**0.375d0
      elseif(INDX.eq.3)then
      rho=1-0.1d0*x2*(1-x2/28.d0*(1-x2/54.d0))
	end if

	return
	end function

c***********************************************************************************************
c///////////////////////////////////////////////////////////////////////////////////////////////

      doubleprecision function find_average(x1,x2) result(average)
c     to find the average value if the integrand 1/2/x^4 in the domain x1-x2
      implicit none
      real(8) x1,x2

      average=0.5d0/(x1*x2)**2

	return
	end function
